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Abstract 

This paper analyzes the nonlinear correspondence between the re- 
flectivity profile (model) and the plane wave impulse response at the 
boundary (data) for a three-dimensional half space consisting of a 
sequence of homogeneous horizontal layers. This correspondence is 
of importance in geophysical imaging, where it has been studied for 
more than half a century from a variety of perspectives. The main 
contribution of the present paper is to derive something new in the 
context of a time-limited deterministic approach: (i) an exact finite 
(non-asymptotic) formula for the data in terms of the model, (ii) a 
corresponding exact inverse algorithm, and (iii) a precise characteri- 
zation of the inherent nonlinearity. Regarding (iii), for generic models 
the correspondence is characterized as a pair of maps, one of which 
is locally linear, and the other of which is locally polynomial. Both 
are determined by a local combinatorial invariant, an integer matrix. 
Concerning (ii), the basic inverse algorithm is modified to allow for 
erroneous amplitude data, taking advantage of the overdeterminacy of 
the inverse problem to recover the exact model even in cases where 
the data is badly distorted. The results are illustrated with numerical 
examples. 
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1 Introduction 

The acoustic reflection problem in the setting of a layered three-dimensional 
half space is fundamental to seismic imaging, a connection in which it has 
been studied for more than half a century. Two recent developments moti- 
vate taking a fresh look at this problem, pared down to its simplest form. 
Firstly, recent progress in superresolution [6] opens up the possibility of 
working with the true impulse response. And secondly, there is an increas- 
ing effort among seismologist to incorporate non-linear effects such as multi- 
ple reflections into the data processing flow [14] [5] [50]. Intuitively, multiple 
reflections are a highly redundant source of potential information about ma- 
terial properties. But the redundancy is wasted without a precise elucidation 
of how material properties are encoded — hence the development of multiple 
suppression (or elimination) as a standard seismic signal processing tech- 
nique [23l Chapter 2]. Indeed, as pointed out the recent survey [15] . seismic 
signal processing is largely predicated on linearization and single scatter- 
ing. The present paper aims to analyze fully the non-linear relationship 
between material properties and boundary measurements for a particular 
formulation of the forward and inverse problems pertinent to a layered half 
space — including a comprehensive treatment of multiple reflections. The 
goal is to establish a mathematically rigorous theory that gives a clear per- 
spective on the deterministic approach. Several new results are established, 
including: 

• exact polynomial formulas for the impulse response (referred to as 
data) in terms of physical parameters (referred to collectively as the 
model) ; 

• the existence of a local combinatorial invariant, an integer matrix, 
that facilitates a precise characterization of the forward and inverse 
mappings between model and data; 

• a fast exact inverse algorithm that avoids downward continuation and 
that exploits the inherent redundancy of the data to correct for erro- 
neous amplitudes. 

The remainder of the introduction is divided into sections as follows. 
Section formulates precisely the forward and inverse problems that are 
the subject of the rest of the paper, and compares the given formulation to 



earlier treatments. Section 1.2 summarizes the paper's main results. And 



Section 11.31 cites some additional related literature. 



12 June 2012 



Peter C. Gibson 



Reflections in layered media 



4 



1.1 Problem statement 

The physical setup is as follows. Consider a three-dimensional acoustic 
medium with coordinates (x, y, z) such that the medium varies only in the 
z-direction and such that the density p{z) and bulk modulus K(z) are piece- 
wise constant with respect to z, having jumps at points 

Z < Zi < ■ ■ ■ < ZM, 

where M > 1. Let Z—\ < zq denote a fixed reference depth in the ho- 
mogeneous half space z < zq from which signals in the form of traveling 
plane waves may be transmitted and received. In keeping with the geo- 
physical perspective, the z-coordinate will be interpreted as depth, and it 
will be depicted as increasing downward, as in Figure [T] The depth range 
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Figure 1: A layered medium, with z increasing downward. 



Zfi—i < z < z n with be referred to as the nth layer, for 1 < n < M. Thus 
there are M layers and M + 1 interfaces, the latter being located at depths 
Zq,Zi, . . . , zm- 

Let u(t, z) denote the velocity (in the z-direction) of a material particle 
at depth z and time t, and let p(t, z) denote the pressure. The medium 
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evolves according to the coupled one-dimensional equations 

du dp . 

1 dp du , 

«i + &= a < llb > 

In analyzing the solution to the above system, there is a choice to be made 
between working with u(t, z) or p(t, z) (or a certain combination or the two). 
The situation is essentially the same whatever the choice; for the sake of 
definiteness the present article will focus on the velocity field u(t, z), which 
is the quantity measured by a coil/magnet geophone, for instance. The 
initial conditions corresponding to a plane wave unit impulse propagating 
downward from z~\ are 

u(0, z) = biz — z-\ ) 

, - (1.2) 

p(0, z) = ^/K[z- X )p{z-x) 5(z - z-i). 

Let G(t) denote the (velocity) impulse response at Z-\, so that 

G(t)=u(t,z-i) (t>0), (1.3) 



the solution at depth Z-i to the system ( l.la|l.lb|l.2 ). 



The work [8j Chapter 3] of Fouque et al. summarizes very clearly the 
standard theory concerning propogation of waves in a piecewise constant lay- 
ered medium, including a derivation of the governing equations ( l.la|l.lb ) 



from physical principles; it will serve as a principal reference. See [SJ Sec- 
tion 2] for an alternate treatment of the same material. The following facts, 
proved in [8j Chapter 3] and elsewhere, serve as a starting point for the 
present paper. For 1 < n < M, let r n denote the two-way travel time (for a 
traveling wave) across the nth layer of the above M-layer medium, and let 
To denote the two-way travel time from depth z-\ to zq. For < n < M, let 
R n denote the reflection coefficient at depth z n relative to a wave traveling 
toward the interface from above. Letting K n and p n denote the density and 
bulk modulus inside the nth layer — with iT_i,p_i and Km+i, Pm+i denot- 
ing the respective values at Z-\ and any point zm+i below zm — the travel 
times and reflectivities are given by the formulas 

^(Zn Z n — 1 

) j d \/K n p n — ^JK n+ ip n+ i 
r n = 7= and Kn = — (1.4) 

yJK n p n V K nPn + V K n+lPn+l 



for < n < M. Note that -1< Rn < 1 by virtue of (jL4j). Let (r, R) denote 
the pair of sequences 

t= (t ,..., t m ) and R = (R Q , R M ), 
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and let |r| denote the total two-way travel time from z-\ to zm, 

\t\ = T H h T M . 

The pair (r, R) will be referred to as a model; it completely determines G(t) 
and hence represents the underlying physical structure insofar as concerns 
the impulse r espo nse. To emphasize this dependence, the impulse response 
as defined in ( 1.3) will be denoted G^ T ' R \t). Moreover, the impulse response 
has the structure of a delta train of the form 

oo 

G ^R) { t ) = Y j a n 5{t-a n ). (1.5) 

n=l 



Without loss of generality the representation (1.5) will be taken to be in 
normal form, whereby each a n ^ and the arrival times a n are in their 
natural order, 

01 < (7 2 < <T 2 < • • ■ . 



In practice the impulse response (1.5) can only be known for some finite 
time interval [0,t max ]. If i max < |r|, say 

T H h TJV < t 

max < Tq + • • • + Tjv+1, 

then the terms rjv + i , . . . , tm and Rn+i , ■ ■ ■ , Rm have no influence on G^ T ' R ^ (t) 
for i G [0, t max ]. That is, letting r' = (ro, . . . , tn) and i?' = (Rq, . . . , -Rat), 

G (r,fl) W=G (r',/?') Wfor t <tmax . 

Thus one may as well study just ,R \ Based on this reasoning, let X[o,|t|] 
denote the characteristic function of the finite interval [0, |r|] and consider 
the partial impulse response 

d 

X[o,M]G (T ' i?) (t) = ^a n <5(t-a n ). (1.6) 

n=l 

Letting (a, a) denote the finite sequences 

cr = (cri, . . . ,cr d ) and a = . . . ,a d ), 
the central problem of the present paper is to analyze the mapping 

{t,R) ^ (a, a). 
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The data (a, a) corresponding to (1.6) is well-defined provided that the 



right-hand side is in normal form, a convention that will be in force from 
now on. 



Forward problem. Given an M-layer model (r, R), for some M > 1, 



determine the data (a, a) for which equation (1.6) holds 



Inverse problem. Given the data (a, a) such that equation (1.6) 
holds for an M-layer model (r, R) , for some M > 1, determine the 
model. 



In the above formulation, the vectors r and R are arbitrary, in the sense 
that any 

reR% +1 and R € (-1, 1) M+1 

comprise a possible model, for any M > 1. Three features distinguish the 
above formulation, setting it apart from earlier treatments of the problem: 

1. the data is a delta train and hence purely singular; 

2. the T n are arbitrary — they need not all be equal; 

3. the data is restricted to finite time. 

The purely singular structure fits the scenario of Candes and Fernandez- 
Granda's work [6] on superresolution. Their results suggest that, within 
certain constraints, it may be possible to extract the function X[0,\t\]G^ t ' R ^ (t) 
from measured data, which one expects to be a convolution of the form 

for some non-sigular wavelet /(£). Deconvolution, which has its own long 
history, will not be considered in the present paper. While the purely singu- 
lar structure is suited to superresolution, it is unsuited to classical methods 
such as that of Gelfand-Levitan, as pointed out by Bube and Burridge [5, 
Section 1.2]. The second item above, that the r n need not be equal, dis- 
tinguishes the present formulation from several influential papers on the 
1-D reflection problem, including that of Berryman-Goupillaud- Waters 
Kunetz [11] and Bube-Burridge [5]. Each of the latter assumes constant 
travel times across layers, which turns out to change the problem substan- 
tially from the generic case in which travel times across layers are unrelated 
to one another. (See below for more on this point.) Lastly, the finite time 
restriction of item (3) makes frequency domain methods unsuitable if one 
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is interested in exact formulas rather than series approximations. This mo- 
tivates Browning's work [1] on extrapolation, for example. Indeed series 
approximations underpin a range of approaches to the reflection problem, 
in both 1-D and higher. See the survey [21] for an overview based on the 
Lippmann-Scwinger equation, and [S] for more recent work along the same 
lines. The present paper avoids series approximations altogether. 



1.2 Overview of main results 

The forward and inverse problems stated in Section |1.1| are solved com- 
pletely, facilitating a precise description of the model-data correspondence. 
The picture that emerges will be sketched below, following some remarks on 
the forward problem. An explicit solution to the forward problem is stated 
in Theorem |2.1| and its corollories. The proof of the theorem — deferred to 
Section [6] — entails a comprehesive analysis of scattering sequences, carried 
out by means of combinatorial arguments tailored specifically to the task. 
The resulting formulas have a very nice structure: the amplitudes a n in the 
data turn out to be polynomial functions of the reflectivity R. These polyno- 
mials are homogeneous, have integer coefficients, and are consequently easy 
to code up. The arrival times a n have an even simpler structure. Thus the 
solution to the forward problem leads to a straightforward exact algorithm, 
Algorithm |5.1| 

It is proved in Section [3] that the mapping from models to data is not 
globally injective; distinct models of differing dimensions can produce the 
same data as one another. This ill-posedness is shown to to stem from 
alignment of the travel time vector r = (to, ... , Tm) with the integer lattice 
Z M+1 , such as when the r n are all equal — a basic assumption of numerous 
papers on the 1-D problem! The ill-posedness is remedied by discarding a 
small set of models (having measure zero in each dimension). As shown in 
Corollary |4.3[ the mapping from models to data is injective on the remainder, 
which are referred to as generic models. 

More than this, the set of M-layer generic models is naturally partitioned 
into open sets X^. On each set, the mapping (r,R) h-» (a, a) is given by a 
fixed formula, which has the following form. Viewing r and a as row vectors, 



there is an integer matrix (Definition 4.3) such that 



a = tA^ (1.7) 

(valid throughout X^). And each amplitude a n has the form p^^ n {R) (also 
valid throughout X^), where p^, n is an explicitly given multivariate poly- 
nomial determined by ip. Thus within each X^ the mapping decouples into 
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two independant maps r i— > a and R i— > a. The integer matrix and 
the polynomials p^ >n change from one open set to another. In fact, the 
integer matrix is an invariant that uniquely labels and determines 
the associated polynomials p^^ n - The mapping from generic models to data 
is a diffeomorphism on each X^, and it is globally injective when X^ ranges 
over all possible M. See Theorem |4.1| 

Concerning the inverse problem, the discrete invariant A^ can be com- 
puted directly from the arrival time data a, where (a, a) is the data corre- 
sponding to some (r, R) £ X^ . In practice the travel time vector r and the 



matrix Aw, are simultaneously extracted from a — see Algorithm 5.2 The 
matrix A^ then serves as a pointer, indicating which formulas Pip^ n {R) are 
associated with which amplitudes a n . This allows one to pick out the am- 
plitudes of primary reflections, which have the simplest formulas, and from 
which the reflection coefficients R n may be rapidly computed. Crucially, it 
also allows one to cross-check the reflection coefficients using the amplitudes 
of multiple reflections (a hugely redundant source) making it possible to 
correct amplitudes that may have been distorted. The idea is implemented 



in Algorithm 5.3 An important feature of the inverse algorithm (Algo- 
rithm 5.2), which depends on the irregular structure of generic models, is 
that it sequentially solves for r from a and then for R from a. (With one 
extra piece of data it is possible to solve the inverse problem for non-generic 
data, but the algorithm does not decouple as in the generic case and is 
therefore much slower.) 

This is not to say that the inverse problem can be solved in a practical 
way without limitation. On the contrary, there is a simple geometric way 



to interpret the formula (1.7) as a projection of lattice points that shows 
that for any precision e, there are orientations of the vector r (and hence 
associated ip) for which successive entries in the arrival time data a will be 
separated by less than e, rendering practical computation impossible. See 
Section [~~ 



The key results in the present paper are: the explicit form of the ampli- 



tude polynomials py, jTl (Theorem 2.1), the decoupled inverse algorithm (Al- 



gorithms 5.2 and 5.3), and the characterization of the mapping from models 



to data as locally linear /polynomial and globally injective (Theorems 4.1 



and 4.2). 



1.3 Related literature 

Newton's paper [15] surveys literature on the reflection problem up to 1980. 
Apart from the papers cited already in previous sections, notable results 
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since 1980 include those of Santosa-Schwetlick [16] and Santosa-Symes [17) . 
The past two decades have seen considerable progress in the study of scaling 
limits for randomly layered media, and time reversal methods. The book 
[8] by Fouque et al. systematically covers the latter topics, which may 
be broadly characterized as stochastic. By contrast, the present paper, 
while starting with the same layered medium as in [8], takes an entirely 
deterministic approach — and the analysis is entirely in the time domain. 



2 Solution to the forward problem 

The main result of this section is Theorem |2.1[ the proof, which entails a 
comprehensive analysis of scattering sequences, is deferred to Section [6j The 
theorem requires a number of technical definitions that will also play a role 
in later sections. 



2.1 Key definitions: lattice set and amplitude polynomial 

The non-negative integer lattice plays a central role in the analysis of 

an M-layer model (r, R); fix notation as follows. A point k G J_ M+l will be 
indexed starting with 0, 

k = (k ,h, k M ), 

consistent with the convention for models. For lattice points k,k' G Z +1 , 
the notation k < k' refers to the entrywise comparison of vectors, 

k < k if and only if k n < k n for each < n < M. 

(This will be referred to as the natural partial order on Z M+1 .) By the same 
token, the minimum of a pair of lattice points is to be interpreted entrywise, 

min{ £;,£;'} = {mm{k ,k' }, . . . , min{k M , k' M }}. 

The symbol 1 G indicates the vector of ones, 

1 = (1,1,...,1), 

and for < n < M, the symbol k n denotes the vector consisting of n + 1 
ones followed by M — n zeros, 

f 1 if < j < n 
K i ~ \ ifn + 1 < j <M • 1 > 

So in particular k M = 1 and A; = (1, 0, . . . , 0); the lattice points k°, . . . , k M 
are called primary vectors. 
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Definition 2.1 (Lattice set) The set of M -layer transit count vectors is 
the subset 

C <£> 

consisting of all k such that: (i) k$ = 1; and (ii) for all 1 < n < M , if 
k n > then k n -\ > 0. Given an M -layer travel time sequence r £ R^ +1 , 
define its lattice set as 

Z T M = {ke£ M \(k,T) < (l,r)}. 

(The term "transit count vector" is a reference to scattering sequences; see 
Section [6j) Note that the lattice set £ T M of a travel time sequence is a finite 
set of lattice points, containing each of the primary vectors k°, . . . , k M . The 
notion of lattice set induces an equivalence relation on travel time sequences: 
declare M-layer travel time sequences r and r' equivalent if 27 M = £ T M . The 
equivalence relation for M = 3 is depicted in Figure [3] of Section |3.2| 

The following (standard) multi-index notation helps to simplify some 
needed formulas. Give a vector x = (xo, • • • , xm ) and a lattice point k £ 
Z M+1 , x k denotes the monomial 

M 

x k = n^-- 

71=0 

And given lattice points b,k £ with < b < k, the multi-index 

binomial coefficient (^) is defined by the formula 

k\ T r I k r/ \ i r ^r?. • 



n r =n 



b J L n\b n J LL n b n \{k n -b n )V 

/ n=0 x 7 n=0 

Definition 2.2 (Amplitude polynomial) Given M > 1, /etx = (xo, • • • ,xm) 
denote a vector of variables; define auxiliary variables y = (yo, . . . ,y&,i) by 
the formula 

y n = y/l-x* (0 < n < M). 
For each k £ £m, let k denote its left shift, k = (fei, k2, ■ ■ ■ , kM, 0). Set 

u = min{l, k} and V(k) = {& £ Z M+1 | u < & < min{/c, k}\ . 

The k-amplitude polynomial a(x, k) is defined as 

ftev(fc) v 7 v 7 
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The following propositions record properties and examples of ampli- 
tude polynomials that will be needed later. Observe that a(x, k) has de- 
gree 2\k\ — 1 and is homogenous with respect to the 2M + 1 variables 
Xo, . . . , xmi 2/0) • • • > VM-l- (Since b < k the variable yu does not occur.) 
Homogeneity and odd degree imply that amplitude polynomials are odd, as 
follows. 

Proposition 2.1 For every k S &m> cl(—x, k) = —a(x, k). 



Proof. Since each y n has even degree in ( |2.2[ ), replacing x with — x is 
equivalent to replacing (x, y) with {—x, —y). Homogeneity implies that this 
is equivalent to multiplication by (— l) de s a ( x ' fe ) = —1. | 

The following connection between amplitude polynomials of differing 
numbers of variables is immediate from the definition. 

Proposition 2.2 LetN>l. Setx = (xo, ■ ■ ■ ,xm) andx' = (xq, . . . ,xm+n)- 
Let k £ £jv/ and k' G £m+n be such that 



Then a(x, k) = a(x' , k') . 



k n ifO<n<M 

ifM + l<n<M + N ' 



Thus an amplitude polynomial a(x, k) depends only on the non-zero entries 
k n of k and the corresponding variables x n . An important special case is 
that of primary vectors k n , as defined above by (2.1). 

Proposition 2.3 Let 1 < n < JA, and let k n £ Z M+1 denote the nth 
primary vector. Then V(k n ) = {/c n } and 

a(x, k n ) = x n (l - xl){l -x\)---{l- x 2 n _ x ). 

A second example serves to illustrate how amplitude polynomials relate 
to one another when their respective transit count vectors differ in a single 
entry. Fix notation as follows: for < n < M, let e n denote the standard 
basis vector in ~Z M+1 , whose entries are 



e n 



1 if j = n 
if j / n 
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Proposition 2.4 Let 1 < n < M — 1, and let k, k' G £m satisfy the condi- 
tions 

k n ^i = k n = k n+ \ = 1 and k' = k + e n . (2-3) 

Then 

a(x,k') _ 
a{x, k) 

Proof. Writing u = min{l,&:} and u' = min{l,A/}, observe that u = u' 
and furthermore that V(k) = V(k'). Also, for each b £ V(k), 

k '\- 9 f k \ fk'-u\_fk-u\ 
b)~ W' \b-u)~\b-u) , 

(-xf- b = (-x n - 1 ){-x)' k - b , and x k '- b = x n x k - b . 
In light of (|2.2|) the conclusion of the proposition follows immediately. I 



2.2 Formula for the data in terms of the model 

The impulse response may be expressed in terms of the underlying model 



using the amplitude polynomials of Definition 2.2, as follows. 
Theorem 2.1 Let (t,R) be an arbitrary M-layer model. Then 

G {r,R) {t)= £ a(R,k)S(t-(k,r)). (2.4) 
fce£ M 

The proof of the theorem is postponed to Section |6j The simplicity of this 



result is remarkable — yet it is apparently new. Restricting the sum (2.4) to 
the lattice set for r adapts the result to the finite time interval < t < |r|, 
as follows. 

Corollary 2.5 Let (t,R) be an arbitrary M-layer model. Then 

X[0M G^ R \t)= Y, a(R,k)5(t-(k,r)). (2.5) 

Thus the data (a, a) corresponding to a given model (r, R) is obtained by 



putting the right-hand side of (2.5) in normal form, thereby solving the 
forward problem. It will be useful to formalize this procedure by introducing 
an explicit mapping, in order to distinguish some fundamentally different 
cases. 
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Definition 2.3 (Enumeration function) Given an M -layer model (t, R), 
set S = suppx[o,\r\]G^ T ' R \t) and d = #S. The enumeration function of 
(t,R), denoted ip = \I/(t, R), is the mapping 

^:Q T M ^{0,l,...,d} 

defined as follows. If (k,r) S, then set ip(k) = 0; otherwise set 

*l>{k) = l + #{ f j G S\a < (k,r)} . 

The enumeration function = ^(r, R) need not be injective; in general it 
induces a partial ordering ■< on the lattice set £ T M according to the prescrip- 
tion 

Vfc, k! G £ T M , k<k' ip(k) < ip{k'). (2.6) 

The data for a given model may be expressed in terms of the model with 
the aid of its enumeration function as follows. 

Corollary 2.6 Given an M-layer model (t,R), write 

d 

n=l 

in normal form. Then, for 1 < n < d, 

a n = (k,r) for any k G ^ _1 (n), 

and 

a n = ^2 a ( R > k )- 

fc6i/' _1 (n) 

This solves the forward problem completely. 

2.3 A geometric interpretation of arrival times 

There is a simple geometric interpretation of the arrival times a n in the 
typical case where ip is non-zero. Up to rescaling, arrival times are essentially 
the orthogonal projection of the positive octant of the M + 1-dimensional 
integer lattice onto a line having direction vector r. More precisely, they are 
a rescaling of the projection of the lattice set £ T M onto the line with direction 
vector r. This can be sketched in low dimensions; as an illustration, consider 
the 3-layer case 

r = (1, 0.327971, 0.152455, 1.51957). 



12 June 2012 



Peter C. Gibson 



Reflections in layered media 



15 




Figure 2: The lattice set £3 (orange lattice points) along with £3 (gray 
and orange lattice points), in the case r = (1,0.327971,0.152455, 1.51957). 
Vectors in £3 have the form k = (1, k\, k%, 0), except (1,1,1,1), which is 
omitted. Points k £ £3 are numbered according to ip(k). 



Except for the vector A; 3 = (1, 1, 1, 1), each member of £3 lies in a two- 
dimensional plane in Z 4 . See Figure [2] For this particular choice of r, the 
numbers a n corresponding to arrival time data are depicted in Figure [2] as 
the orthogonal projection (red dashes) of the orange points onto the blue 
line. The blue line is the orthogonal projection of r onto the (ki, A;2)-plane. 
Later arrival times would correspond to projections of the gray points. The 
numbering of the orange points indicates the values of the enumeration 
function which will be shown later to be a locally constant function of r. 

This view of arrival times as the orthogonal projection of lattice points 
onto a line explains both their pattern and density with increasing time. It 
is closely related to the cut-and-project construction of quasicrystals due to 
Kramer and Neri [IU| . In this connection see [1] and also [7j, who traces 
lattice projections back to the almost periodic functions of Bohr [3]. Note 
that for layered media the link to almost periodicity has a direct manifesta- 
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tion: the Fourier transform G^ t ' R ^(lo) of the impulse response is an almost 
periodic function in the sense of Bohr. Returning to present considerations, 
as long as the direction of the line determined by r is not rational, arrival 
times will be aperiodic and eventually arbitrarily dense. On the other hand, 
the periodic structure of arrival times associated with, for example, a Bragg 
mirror, corresponds to a rational choice of r that lines up with the lattice. 



3 Ill-posedness of the inverse problem 
3.1 Examples of non-injectivity 

The mapping (t, R) \— > (a, a) taking a model to its data is not injective, as 
the following example illustrates. 

Example 3.1 Consider a 2-layer model of the form 

(t,R) = ((t ,ti,ti), (R , -j=,Rq)), 
and the accompanying one-layer model 



(T',R') = [(ro,r 1 ),(Ro,^=)). 



Evaluation of the formula in Corollary \2. 5| shows that 
X[o,M]G (T ^(t) = X[o, H ] G^')(t) = R S{t - r ) + ^ 6[t - (r + n)). 

Thus the two distinct models give rise to the same data. 

Let t/j = *$>(t,R), the enumeration function for the above 2-layer model. 
Notice that 

^((1,2,0)) =i;(k 2 ) =0. 



So ip both takes the value zero and fails to be injective. Example 3.1 is part 
of a general pattern; similar examples exist in every higher dimension. The 
following technical lemma captures the essential reason. 

Lemma 3.2 Let r be an M + 1-layer travel time vector such that the set 
S, consisting of all k G £]vf+i with the property that k ^ k M+l and (k, r) = 
(k M+1 ,r), has at least one element. Let 1Z denote the set of all points 
(Ro, . . . , Rm) £ (— 1, 1) 1 for which there exists a non-zero value Rm+i G 
(— 1, 1) such that 

a{R',k M+1 )+Y,a(R',k) = 0, (3.1) 

kes 
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where R' = (Ro, . . . ,Rm+i)- Then the Lebesgue measure of 1Z in R M+1 is 
strictly positive. 

Proof. Note that if k G is different from fc +1 , then a(R',k) is 

independent of Rm+i- Moreover, letting n denote the maximal index such 



that k n ^ 0, the formula (2.4) shows that a(R',k) is a multiple of R n . 



Therefore the set V of points 

(i?0,... J i?M)G(-lAl/2) M + 1 



< (3/4) 



M+l 



for which 

< ^a(R',k) 

has positive measure, since ^ SfceS a (-^'' ^) ~ ^ as maxo< n <M \Rn\ — > 0. 
It follows from Proposition 2.3 that for each (Rq, . . . , Rm) £ V, setting 

M 

R M+1 = -*£a(R>,k)/ll(l-R 2 n ) 



kes 



n=0 



yields the desired equation (3.1). 



Many of the historical approaches to the 1-D seismic problem make a 
simplifying assumption that the travel time vector r is constant. But con- 
stant t are easily seen to satisfy the hypothesis of Lemma |3.2| It follows 
from this that the inverse problem for such models is ill-posed. 

Theorem 3.1 Let r be a constant M -layer travel time vector, and let r' de- 
notes its extension to an M+l-layer constant vector. Write 1Z C (—1, 1) M+1 
for the set of all points R = (Ro, • • • , Rm) f or which there exists a non-zero 
Rm+i G ( — 1, 1) such that the models (t,R) and (t',R') generate the same 
data, where R' = (Ro, . . . , Rm+i)- The set 1Z has positive Lebesgue measure 
in R M+1 . 



Proof. Observe that r' satisfies the hypothesis of Lemma |3.2[ consider the 
reflectivities R' supplied by the lemma. Note that £^ +1 \ {^ A/+1 } consists 
of the elements of £ M , each extended by a entry; therefore the supports 
of 

X[0,]r|]G^(t) and X[o,\r'\]G^' R 'Ht) 

differ at most by (k M+1 , r'). But the latter does not belong to the suppo rt of 
X{o,\t'\]G( t ' > R >(t) , since its coefficient is given by the left side of equation (3.1 ) 
which, by Lemma 3.2, is zero. The remaining coefncents are identical be- 



tween the two impulse responses, so (r, R) and (r', R') generate the same 
data. I 



12 June 2012 



Peter C. Gibson 



Reflections in layered media 



18 



3.2 The hypothesis of Lemma |3. 2 



What about non-constant travel time vectors? After all, a randomly chosen 
M-layer travel time vector will not be constant. (In fact, with probability 
one, it's entries will be linearly independent over the rational numbers.) 
A starting point is to consider travel time vectors that satisfy the hy- 



pothesis of Lemma 3.2 since this is the source of the ill-posedness expressed 
in Theorem 3.1. The case of 3-layer travel time vectors r = (to, t\, t 2 , 73) 
illustrates the general situation. Observe that a transit count vector k £ £F M 
different from A: 3 satisfies (k, t) = (/c 3 , t) if and only if 

(k% - l)fi + (k 2 - 1)t 2 - t 3 = 0, where f = n+ 4 +75 ( r i, r 2 , r 3 ). (3.2) 
Thus the normalized subvector f determines whether a given 3-layer transit 



vector t satisfies the hypothesis of Lemma 3.2 Of course f is a point in 



the standard two dimensional simplex A 2 , which is a triangle. Each k S £3 



determines a line segment in A 2 defined by the equation (3.2). The set of 
all such lines is plotted in Figure |3| The vertices, clockwise from the top, 
are f = (0, 0, 1), (0, 1, 0) and (1, 0, 0). The union of lines represents the set 



of transit count vectors that satisfy the hypothsis of Lemma 3.2, so the 
lemma fails to hold on the complement. The cells bounded by the lines 
represent the set of all travel time vectors that have a given lattice set, in 
the sense that £3 is constant for f ranging over a single cell — with distinct 
cells corresponding to distinct lattice sets. The bottom left triangular cell 
represents the set of travel time vectors t such that ££ = {k°, k , k 2 , A: 3 }, the 
minimum possible. Moving to the right, each successive cell adds a vector 
of the form k = (1, n, 0, 0), with n increasing from 2 to 00. Moving upward 
from the lower left, each successive cell adds a vector of the form (1, 1, n, 0). 
A zoomed image from the upper part of the triangle is depicted in Figure |4| 
Of course there is a corresponding construction in every higher dimen- 
sion. In general, for M-layer travel time vectors, the triangle A 2 becomes a 
standard [M— l)-dimensional simplex Ajv^—i, lines become hyperplanes, and 



polygonal cells become convex polytopes. Lemma 3.2 fails on the interior of 
each of the convex polytopes. 

Thus in general Lemma |3.2| only holds on a set of M-layer travel time 
vectors of measure zero in K >0 + ; this set happens to include constant travel 
time vectors. The idea of the next section is to recover injectivity of the 
mapping from models to data by excluding a thin set of models which 
encompasses all those whose travel time vectors satisfy the hypothesis of 
Lemma 13.21 
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Figure 3: The 3-layer travel time vectors that satisfy 



Lemma 3.2, represented by line segements within the tri- 



angle. Each cell bounded by line segments represents the 
set of all vectors having a common lattice set. 



4 The inverse problem for generic models 

4.1 Definition and characterization of generic models 

Definition 4.1 (Generic) A model (t,R) is termed generic if its enumer- 
ation function is injective and never zero. 

Thus if (r, R) is generic, then its enumeration function ip = ^(r, R) has an 
inverse, 

i P = ^- 1 :{l,...,d}^Z r M . (4.1) 
Here is an alternate characterization that follows directly from Defintion|2.3| 



Proposition 4.1 An M-layer model (t,R) is generic if and only if the 
following conditions hold. 
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Figure 4: A zoomed portion of the upper part of the 
triangle in Figure [3} The green + indicates the location 
of the travel time sequence depicted in Figure [2] 



1. The mapping k i— > (k,r) is injective on 27 M . 

2. For each k £ £ T M , a(R, k) ^ 0. 

It will be useful to refer to the reflectivity sequences that satisfy property 
(2) of Proposition 4.1 and furthermore to have a notation for the set of 
travel time vectors associated to a common enumeration function. 

Notation 4.2 Let if) = ^?(t, R) for some generic M -layer model (r, R). Set 

U^, = {R' E (-1, 1) M+1 | a(R', k)^0\/ke £ T M } . 
In addition, set 

= |V G Ri J +1 | V R! € #(t', R') = V} • 

Observe that t' S only if both = £7 M and condition (1) of Proposi- 
tion 4.1 holds. More than this, 

Proposition 4.2 If ip = ^(t,R), then LO, C R M+1 is a convex open set 
containing t, and IZ^p C (—1, is an open set containing R. 
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As in Section 3.2, the case M = 3 serves to illustrate the general situation. 



Proposition 4.1 implies that the question of whether or not a given 

r = (t ,...,t 3 ) G R* 
is generic depends only on the normalized vector 

f = T 1+ r 2 +r 3 ( T l> T 2> r 3). 

In detail, r fails to be generic if and only if there exist k, k' G £3 such that 

((ki,k 2 ,k 3 ),r) <(l,r), ((k[,k' 2 ,k' 3 ),r) < (1,t) (4.2a) 

and (Oi - k[,k 2 - k' 2 , k 3 - k' 3 ),f) = 0. (4.2b) 



Note that the two inequalities (4.2a) simply express that k,k' G £3. The 



particular case where k G £3 and k' = /c 3 reduces to equation (3.2). There- 
fore the set of f that satisfy ( |4.2[ ) for some k, k! G £3 is a superset of the line 
segments plotted in Figure |3j The full set of solutions is plotted in Figure [5j 
with the additional segments coloured orange. The sets are the travel 
time vectors r that correspond to the interiors of the cells in this figure. A 
zoomed image from the upper part of the triangle is depicted in Figure [6j 

A similar picture exists in every higher dimension. Thus in general, 
is the interior of the postive cone over a convex polytope crossed with K>o- 
The important point is that is a convex open neighbourhood of r, and 
such neighbourhoods partition the set of all travel time vectors belonging to 
generic models. 



As for reflectivities, given a fixed ift = ^(t,R), Proposition 4.1 asserts 
that the vectors R' G (— 1, 1) M+1 excluded from are precisely those 
satisfying a polynomial equation of the form 

a(R\ k) = for some k G £ T M . 

Thus TZ^ is the complement of a finite union of algebraic hypersurfaces in 
(-1, 1) M+1 . In particular, is open and has full measure. 

The enumeration function for a generic model may be represented as a 
matrix, as follows. 

Definition 4.3 (Enumeration matrix) Let (r, R) be a generic model, so 
that t/j = \I/ (r, R) is a bijection of the form 
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Figure 5: The 3- layer travel time vectors that fail to 
be generic, represented by line segements within the tri- 
angle. The sets U$ are represented by the open cells 
bounded by line segments. 



The enumeration matrix representing tp is defined to be 

^ = (r 1 (i) T r l vf ■■■ r\df). 

Thus Ajjj is the (M + l) x d integer matrix whose nth column is the transpose 
of the vector ^~ 1 (n) G 



Note that since each of the primary vectors A: , 
matrix 



k M belongs to £? M , the 



K 



(ii i ... i\ 

1 1 ••• 1 
1 ••• 1 



\ 



(4.3) 



1 / 
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Figure 6: A zoomed portion of the upper part of the 
triangle in Figure [5] Again, the green + indicates the 
travel time sequence of Figure [2] 



is a submatrix of A^. Therefore A^ has rank M + 1. For future reference, 
let Jm denote the inverse of Km, 



( 1 
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(4.4) 



4.2 The model-data correspondence 

The notation developed in the previous section sets the stage for a very 
simple description of the correspondence between a generic model and its 
data. 

Theorem 4.1 Fix an enumeration function ip, and let (t,R) G x 7Z^ 
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have associated data 

(a, a) = ((cri, . . . ,a d ), (ori, . . .,a d ))- 
Then the data is given in term of the model by the formulas 

a = tA^ and a n = a(R, (1 < n < d). (4-5) 



The crucial point is that the formulas ( |4.5[ ) hold locally; the same formula is 
valid throughout the open set x 1Z^. And open sets of the form x TZ^ 
partition the set of all generic M-layer models. The inverse mapping, as it 
applies to data for generic models, may be characterized explicitly as follows. 

Theorem 4.2 Let (a, a) be the data corresponding to an M-layer model 

(t,R) e^x K^, 

for some enumeration function ip, and let a^p denote the subvector 

o~i\> = (o>(fc°)' ^(fe 1 )' • ■ ■ >°>(fc M )) G ^ +1 - 
Then the model is given in terms of the data by the formulas 

t = o-^Jm, (4.6) 
Ro = ai , Rn = n) (l<n<M), (4.7) 

11.7=0 1 1 n j) 



where Jm is the matrix (4-4), and, as usual, k n denotes the nth primary 
travel time vector. 

Thus the restriction of inverse map to the image of x 1Z^ is linear in the 
arrival time data, by (4.6). The recursive scheme (4.7) captures the precise 
non-linearity of the reflectivities, as determined by the amplitude data. An 
important consequence of the inverse formulas is injectivity of the forward 
map. 

Corollary 4.3 Restricted to generic models, the mapping from models to 
data is injective. 



Proof Local injectivity of the forward map is implicit in Theorem 4.2 



which gives the pre- image of a point in the range of a given patch U$ x 7Z^. 
Global injectivity of the mapping (r, R) \— > (<r, a) follows from the inverse 



12 June 2012 



Peter C. Gibson 



Reflections in layered media 



25 



algorithm in Section [5j Algorithm 5.2, which computes r — and hence ip- 



directly from a. Thus if two models have the same image, they must belong 
to the same patch and therefore be identical. I 



By Corolloary 4.3 the question of well-posedness of the inverse problem 
boils down to continuous dependence of the model on the data. If one re- 
stricts the question to the image of a particular open set of generic models 
U-ij) x TZ^p, there is no problem. The forward and inverse maps are diffeomor- 
phisms. Considered globally, however, the situation is much more delicate. 
To take a concrete example, consider a point 

a = (<ri,f7 2; 0-3,0-4) G K 4 

as a candidate for travel time data. Is a in the range of the forward map- 
ping? And if so, what is the structure of the set of nearby data points? By 



Theorem 4.2 , the set of arrival time vectors a £ R d belonging to the range of 
the forward map is a finite union of convex polytopes of various dimensions 
< d, each of which is open relative to its affine hull. The case d = 4 is 
illustrated in Figure [7j where the data is projected onto the (02, <73)-plane, 
and normalized so that it lies inside the unit square. There are two separate 
pieces, a triangle, which corresponds to 3-layer models whose data has 4 
terms, and a segment, which corresponds to 2-layer models whose data has 
4 terms. The intersection of the closures of these two pieces is a single point 



that is just like Example 3.1 It is clear even in this low dimensional case 
that an arbitrary point a G R 4 will not in general have a unique closest 
point in the closure of the range of the forward mapping. Thus there is no 
clear-cut way to project a given candidate for a data point onto the range. 



(See Section 5.2.1 for a different approach.) Thus the inverse problem is 
well-posed if one restricts to the range of the forward mapping; but the 
range itself has an irregular structure. 



5 Algorithms and examples 

With an explicit formula for the truncated impulse response 

d 

X[o,M]G (T ' R) (t) = ^a„5(t-a n ) 

n=l 

in terms of the underlying model (r, R) , it is a straighforward matter to com- 
pute the exact impulse response efficiently, provided the number of layers is 
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0.4 - 



0.2 - 



0.2 0.4 0.6 0.8 1.0 



Figure 7: The set of arrival time vectors a € R 4 belonging 
to the range of the forward map, normalized and projected 
onto the (o"2, o"3)-plane. The triangle is the image of 3-layer 
models, and the segment is the image of 2-layer models. 



not so large as to render the inherent polynomial evaluations prohibitively 
expensive. The corresponding inverse algorithm, while perhaps less obvious, 
does not require the polynomial evaluations inherent in the forward algo- 
rithm. Rather, it is constrained by search and sort procedures, details of 



which are given in Section 5.2 



5.1 Forward algorithm 

The following algorithm encodes the solution to the forward problem, using 



formula (6.14) of Theorem 2.1 



Algorithm 5.1 (Forward) 

Input: (t, R) = ((t q , . . . , t m ), (Rq, Rm)) 
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Step 1: Allocate £ T M = (jfeW, . . . , fcW) . 

Siep 2: Evaluate r ■ £ T M = ((r, fcW), . . . , (r, . 

Step 3: Sort r ■ S7 M into increasing order, a = f(r, k^ 1 "), . . . ,(t, fc^^}) 
Step 4'- Compute a n = a(R, k^( n '') for n = 1, . . . ,d, using formula (6. 14) 
Output: (a, a) = ((at, . . .,a d ), (a%, . . .,a d )) 

The output of the algorithm is the arrival times and amplitudes for the 
normal form of G^ T ' R \t), restricted to t < |r|, 



X[o,[r[]G (T ' fl) (*) = J2a n S(t-a n ). 



n=l 



5.1.1 An example with 17 reflectors 

Figure M depicts a 17-layer model (where R n is plotted against tq + • • • 



r n ), and its impulse response computed by Algorithm 5.1 (in approximately 



17.1 seconds). This impulse response has 33714 terms, most of which are 
low amplitude late arrivals, invisible in Figure § (b). Note the very small 
amplitude 0.012 term at 8.75 seconds of the model. While the first four 
terms in the model have easily recognized corresponding primary reflections 
in the impulse response, this small amplitude is totally obscured. There is 
no way to recognize its associated primary reflection by simply studying the 
plot of the impulse response, no matter what the scale. 




I i 



A a, 



Figure 8: A 16 layer model (i.e. with 17 reflectors), (a) The true model, (b) 
The corresponding impulse response computed by Algorithm 5.1 



Figure [9] (a) shows a zoomed view of the impulse response (plotted with- 
out stems) revealing the explosion of later arrivals. Figure^] (b) shows a 7003 
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term decimated version of the impulse response in which all terms having 
absoute amplitude < 10~ 4 have been removed. 
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Figure 9: A zoomed version of the impulse response, (a) The true impulse 
response (plotted without stems), and (b) a decimated version. 



5.2 Inverse algorithm 

For the inverse algorithm, a further subdivision of £ T M is needed. Given 
T (iy — ( T0) ^ Tn ~j an j s > o, write 

£;J 5 = |fe e £ n fc n > 1 and ( T W, fc) < s j . 

Algorithm 5.2 (Inverse) 

Jnpui: (<r, a) = ((cri, . . . , a d ), (ai, . . . , a d )) 
Stage I: Arrival time inversion 

Step 1: Set (to, ti) = (o"i, 02 — 0i), S 1 = {ex.,- | 3 < j < d} and n = 1. 

5tep £a: Given = (tq, . . . ,T n ), allocate ££^.. 

Step 2b: Set S = S \ { (r^,k) \ k e tfS = ®, go to Step 3. 

Step 2c: Set T n+ \ = minS" — \t^\ and n = n + 1; return to Step 2a. 
Step 3: SetM = n, r = tK 
S^age 77: Amplitude inversion 

Step 4-' For each < N < M , determine p{N) such that cr p (7v) = (r, ). 

Step 6: For each 1< N < M, set R N = — ° p(JV) ;f N ~t , . 
Output: (r, 7?) = ((r , . . . , r M ), (i?o, • • • , -Rm)) 
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Note that Step 6 uses the result detailed in Proposition |2.3[ The above 



inverse algorithm has a number of remarkable properties, as follows. 

• The algorithm is exact. 

• The arrival time inversion is independent of amplitude data. 

• If the input arrival times are collectively shifted by k, giving modified 
input 

(a + re, a) = ((<ri + k, . . . , a d + re), (ai, a d )) , 
then the resulting output is (f, R) = ((to + k,ti, . . . , tm), (Rq, ■ ■ ■ , Rm))- 

• Given any subsequence of the original input data that contains all the 
primary reflections, 

(a, a) = ((er mi , . . . , cr m J, (a mi , . . . ,a m J) , 

the resulting output is (r,R). In other words, the algorithm recovers 
the model exactly given only partial data. 

• The algorithm is fast, and can be accelerated by deleting low amplitude 
data, provided no primary reflections are deleted. 

The inverse algorithm recovers the 16-layer model depicted in Figure [8] 
from its impulse response exactly; the computation took 127 seconds in 
a less-than-optimal implementation in Mathematica. Remarkably, Algo- 



rithm 5.2 also correctly recovers the model from its decimated impulse re- 
sponse depicted in Figure |9j in this case the computation took just 28.2 
seconds. In both cases the recovered model is identical with the plot in 
Figure |8^a) (including the small amplitude term at 8.75 seconds!). 

Figure [10] depicts the impulse response for a 7-layer model, with some 



randomly generated spurious arrivals added on top. While Algorithm 5.2 
recovers the model exactly from the uncorrupted data, it falters when fed 
the data that includes spurious arrivals. See Figure |TTj 

5.2.1 Modification to detect spurious arrivals 

The inverse algorithm can be modified to handle the presence of spurious 
arrivals. One way to control for this is the following modification to Step 2a. 



Step 2a(i): Given rW = (tq, . • • , r n ), allocate £ 



_(n) 

n,a d m 
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(a) 



(b) , 



Figure 10: (a) The true impulse response for a 7- layer model (115 data 
points), and (b) with 12 spurious arrivals added. 



(a) 




(b) 



' f ' i ■ ■' i 



6 8 



Figure 11: (a) The true 7-layer model corresponding to Figure 10 (b) The 
partial model computed from by Algorithm |5.2| with the corrupted data as 
input; the computation was terminated after 612 seconds. 



Step 2a(ii): If S n { (r^ n \k) \ k E = {minS} / S, then 

set S = S \ {mills'}, set r n = minS 1 — |r^ n_1 ^|, and return to Step 2a(i 

The effect of this modification is to delete data points that appear to cor- 
respond to primary reflections, but for which there are no corresponding 



multiple reflections in the later data points. The modified Algorithm 5.2 



with the corrupted version of the above impulse response as input, recov- 
ers the model in Figure [lT|(a) exactly. The computation took just 0.0569 
seconds. 



5.2.2 Modification to correct for distorted amplitudes 

Further, one can exploit redundancy of multiples to control for inaccurate 



amplitude data. The idea is to use the formula (2.2 ) to solve for products and 
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ratios of reflection coefficients that are given by distinct data points. There 
are many ways to do this. As a simple illustration of the idea, the following 
scheme uses Proposition |2.4| Suppose that the above inverse algorithm has 
input (a, a), where the amplitude a is distorted, and let (r,R') denote the 
resulting output — where the ' indicates that the reflectivity is distorted. 
Assuming that r is accurate, the following third stage is designed to follow 



Algorithm 5.2 



Algorithm 5.3 (Reflectivity correction) 

Input: (r, R') = ((to, . . . , tm), (R'o, • • • , R'm))> an< ^ (°"> a )> the respective 

output from, and input to, Algorithm 5.2 
Stage III: Reflectivity correction 

Step 7: For each 1 < n < M — 3, allocate the set E n of all pairs (k, k') 



from £7~ M satisfying condition (2.3) in Proposition 



2., 



Ct I 

Step 8: For each 1 < n < M — 3, compute the set C n of all ratios 
such that o-j = (k,r) and a.y = {k',r), where (k,k') £ E n . 

Step 9: For each 1 < n < M — 3, set c n to be the (mean of the) most 
common value(s) in C n . 

Step 10: Set R' ' = R' . For n from 1 to M — 4, set R^ = CnjR" n _ x . 

Step 11: For n from M — 3 to M, let j n be such that aj n = {k n ,r) (with 



k n as in Proposition 2.3), and set R" n = a jn /Y[^(l - (R") 2 ) ■ 
Output: R" = (Rl,...,R" M ) 

The above reflectivity correction requires that the input value R' Q be 
accurate — meaning that R' = Ro, the true reflectivity — and that the later 
amplitudes a Jn invoked in Step 11 also be accurate. As long as these assump- 



tions are satisfied Algorithm 5.3 can recover the true reflectivities. Figure 12 
shows an example. The short-term sine wave distortion of the amplitude is 
meant to mimic the phenomenon of ground roll endemic in land-based seis- 
mic acquisition. 



The results of Algorithm 5.2 followed by Algorithm 5.3 are depicted in 
Figure 13 Despite the substantial sine wave distortion, the exact reflection 
coefficients are extracted using the set C n . Without Algorithm |5.3| only the 
first four reflection coefficients would have been correctly estimated. 
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Figure 12: An 11 layer model with distorted amplitudes, (a) The true im- 
pulse response (27108 data points), (b) A 0.2 amplitude sine wave distortion 
added to the true impulse response between t = 2.51s and t = 6.56s. The 
sets Ci , . . . , Cg are plotted for the true impulse response (c) and the distorted 
impulse response (d). The range of distortion is clearly evident. 



6 Proof of Theorem 2.1 



A standard idea is to view G<- T < R \t) as the sum total of all possible sequences 
of successive reflections and transmissions of the initial pulse that eventually 
return to z~\. The sequence of transmission coefficients T = (To, . . . ,Tm) 
corresponding to a given sequence of reflection coefficients R = (Rq, . . . , Rm) 
is given by the formulas 

T n = y/l-Rl (0 < n < M). 

(See P, Chapter 3] for details.) The idea is worth illustrating, since it 
underpins the arguments below. For example, after tq/2 seconds, the initial 



downward traveling pulse 5(z — Z-\ — Cq€) satisfying equations ( l.la|l.lb|1.2 ) 
reaches the interface at zq. Part of it is transmitted into the first layer as 
Tq5(z — zq — c\{t — y))' which, after another n/2 seconds, reaches the 
interface z\. Part of this pulse is reflected back into first layer as RxT^biz — 
z\ + c\(t — T °^ T1 )) • Traveling back up to Zq, the latter reaches zq after a 
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Figure 13: Amplitude correction applied to the impulse response for the 11 



layer model of Figure 12 The true model (a) and the model computed by 



Algorithm |5.2| (b). The correction after Step 10 of Algorithm 5.3 (c) and 
after Step 11 (d). The true model is recovered exactly. 



further n/2 seconds, and is partly transmitted back to Z-\ as 
R^diz-zo + coit-^^)), 



(6.1) 



arriving at z~\ at time a = tq + r\. Thus the part of the initial pulse 
that traverses the sequence of depths p = (z_i, zo, z±, zo, z_i) returns to 
Z-\ with an amplitude a = R\Tq at arrival time a = tq + t%, thereby 
contributing a term of the form a5(t — a) to G^ T ' R \t). The sequence p is 
called a scattering sequence, and the amplitude a is called the weight of p. 
The impulse response itself is a delta train of the form 



(6.2) 



n=l 



composed of the cumulative contributions of all possible scattering sequences 
returning to with their associated weights and arrival times. In the 
above example, p is the only scattering sequence having arrival time a. But 
in general different scattering sequences may arrive simultaneously, so each 
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amplitude a n occurring in (6.2) is the sum of the weights of all scattering 
sequences arriving at time a n . 

While the interpretation of the impulse response in terms of scattering 
sequences as outlined above is well established, it has not been used as a 
means to express G( T,R \t) directly in terms of the model (t,R). This is 
perhaps because a full analysis of all possible scattering sequences is viewed 
as being enormously complicated, with any resulting representation being 
hopelessly cumbersome and of little practical use. Such a view seems unjus- 
tified. It is shown in the present section, that, with the right perspective, 
a direct analysis of scattering sequences leads to t he surprisingly simple, 
explicit formula for G^ T ' R \t) stated in Theorem 



2.1 



6.1 Scattering sequences and the impulse response 

The first step is to establish some terminology and notation, as follows. A 
scattering sequence that starts and ends at Z-\ is represented by a path in 
the graph 



Z-i Zo Z\ Zi Zi Za --■ Zm-] Zm 

(6.3) 

Formally, given an integer M > 1, let Sm denote the set of all sequences of 
the form 

P = (Po, Pi, • • • , Pi) 

such that L > 2 and: 

Po = PL = Z-\ and p„ £ {z\, . . . , zm} if 1 < n < L — 1; (6.4a) 

Vn with < n < L — 1, 3j with — 1 < j < M — 1 such that 

" " J ~ J - (6.4b) 

|P„, p„+lj = {Zj,Z j+ l\. 

The elements of Sm will be referred to as scattering sequences. Condition 



(6.4a) says a scattering sequence starts and ends at z~i, and condition (6.4b) 



(which refers to unordered pairs) says that adjacent terms in a scattering 



sequence are adjacent vertices in the graph (6.3). 

For example, for < n < M, the shortest scattering sequence that 
reaches z n is 

(z-i, z ,zi,..., z n -i, z n , z n -i, ... ,zi,z , z-i); (6.5) 

this is called a primary scattering sequence. A scattering sequence reaching 
maximum depth z n that is not shortest possible is called a multiple scattering 
sequence. 
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6.1.1 The weight of a scattering sequence 

The weight w(p) corresponding to a scattering sequence 

P = (Po,Pi,---,Pl) 

in an M-layer model (r, R) is defined as follows. For each n in the range 
1 < n < L — 1, and given that p n = Zj, define 

Rj if p n _l = p„+l = Zj-l 

-Rj if Pn-l = Pn+i = Zj+l . (6.6) 



R?- otherwise 



The three possibilities correspond respectively to: reflection inside the jth 
layer at Zj; reflection inside the (j+l)st layer at Zj\ and transmission between 
the jth and (j + l)st layers. Finally, set 



L-l 



^(p) = n Wn - 



n=l 

Thus the part of an initial unit impulse that traverses p returns to Z-\ with 
amplitude w(p). 

6.2 Transit count and branch count vectors 

The purpose of this section is to define two maps, 

k,/3:S m ^Z M+1 

that associate integer vectors to a given scattering sequence. 

A scattering sequence p G Sm rnay be represented graphically as in 



Figure 14; Stanley [18] calls such a representation a Dyck path. 

Given the Dyck path for a scattering sequence p 6 Sm, let t denote the 
horizontal coordinate and (as usual) let z denote the vertical coordinate. 
Let U denote the portion of the t, z-plane on or above the Dyck path and at 
or below z = Z-\ — the shaded region in Figure fT5| For each n in the range 



< n < M, consider the horizontal line L n at depth z n . The intersection 
L n DU consists of a disjoint union of closed intervals J™, where 1 < j < k n ; 



see Figure 15 The intervals of If are of two types: degenerate intervals 



consisting of a single points; and non-degenerate intervals having positive 
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length. Letting k n denote the total number of intervals and letting b n 
denote the number of non-degenerate intervals, set 

«(p) = k = (k , . . . , k M ) and /3(p) = b = (b , . . . , b M )- 

Observe that the entry k n of the vector k = k(p) counts the number of 
times the Dyck path crosses back and forth across the nth layer z n _\ < z < 
z n ; the vector k is therefore called the transit count vector for p. The vector 
b = /3(p) is called the branch count vector, for reasons that will be apparent 
in Section lOl 

6.2.1 The set of all transit count vectors 

The range of k is contained in the set 

-Cm = 

{(fc , h,..., k M ) e N M+1 | A; = 1 k Vn < M - 1, k n = => k n+1 = } . 

(6.7) 

Conversely, for any given k £ £m, it is straightforward to construct a re- 
alizing scattering sequence. Thus £m is precisely the range k(Sm) of the 
mapping k : Sm — > Z M+1 . 
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6.2.2 The formula for arrival times 

The arrival time of p, expressed in terms of the transit count vector k = n(p), 
is simply 

(k, t) = k T + km H h k M T M - 

Therefore G^ R \t) may be written as 

G (r,R) {t)= £ a(R,k)S(t-(k,r)), (6.8) 
kez M 

where the amplitudes a(R, k) are given by the formula 

a(R,k)= Y, w b)- ( 6 - 9 ) 

P6Sm such that 
«(p)=fc 

Note that the weight w(p) of a scattering sequence depends only on R, and 
not on t; the next step is to find an explicit formula for a(R,k). This 
is greatly facilitated by introducing another representation for scattering 
sequences, in terms of trees. 

6.3 The tree representation of a scattering sequence 

A tree is a connected cycle-free graph. The vertices of a tree are divided into 
three anatomical types, as follows: the root is a single, specially designated 
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vertex; a non-root that belongs to just one edge is called a leaf; all other 
non-root vertices are called branch points. Vertices in a tree have a height, 
determined by their distance (in the sense of shortest path) to the root. See 



Figure 16 




Figure 16: A tree, with the leaves coloured green. The root is 
at the bottom, and horizontal lines indicate the various heights 
of vertices. 

The association between a scattering sequence and a tree is well-known 
(see [18} Exercise 6.19]) and arises, for instance, in the analysis of Brownian 
excursions and superprocesses (see [131 Section 1.1]). The tree representing 
a scattering sequence p S Sm m &y be obtained simply by collapsing its Dyck 
path, as follows. Recall the intervals Ij 1 used to define /c(p) and /3(p) above; 
for present purposes let I^ 1 denote the intersection of z = z_i with U. To 
collapse the Dyck path, contract each of the intervals J" to a point, keeping 
the distances between intervals unchanged, and interpolate this horizontal 
contraction linearly on each depth interval < z < z n . This operation 
transforms the original Dyck path into a tree (in fact it is an isotopy be- 
tween the region U and the resulting tree). See Figure 17 Note that the 



degenerate intervals I? (coloured green for emphasis) end up as leaves, while 
non-degenerate intervals are contracted to branch points of the tree — except 
for 17 , which is contracted to the root. 

Conversely, given a tree, one may recover the original scattering sequence 
by tracing the outline of the tree, from the root (keeping the tree on the 
left), and recording the depths of the vertices in the order that they are 
passed. 
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Figure 17: The collapsing of a Dyck path (a) to a tree (f) (depicted upside 
down), an operation which is reversible. 



Observe that the vectors (k,b) = («(p),/3(p)) have a simple interpreta- 
tion in terms of the tree representing p £ Sm- For each < n < M, k n is 
the number of vertices at depth z n , and b n is the number of branch points 
at z n . (This is the reason for calling b the branch count vector for p.) There 
are some evident constraints on these quantities. Note first that bu = 0. 
Furthermore, for < n < M — 1, 



min{l, k n+1 } < b n < k n+1 (0 < n < M - 1) 



(6.10) 



since each vertex at z n +i is connected by an edge to a unique branch point 
at z n . It is convenient to refer to the left shift k of k = k(p); that is, 



k=(k 1 ,k 2 ,...,k M ,0)GZ M+1 . 



In terms of this notation, the constraints (6.10) become 



min{l, k n } < b n < min {k n ,~k n } (0<n<M). 



(6.11) 
(6.12) 
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An easy application of the tree representation is to determine the possible 
values of b = f3(p) for each given k = ft(p), or in other words, to determine 
the set 

V(k) = f3(K~ x (k)). (6.13) 



In fact V(k) is determined precisely by (6.12). 
Proposition 6.1 Given k G £m? 

V(k) = {b | min{l, k} < b < min{/c, k}}; 
equivalently, V(k) may be expressed as a Cartesian product of sets, 

V(k) = V xV 1 x---V M , 

where V n = |min{l, k n }, min{l, k n } + 1, . . . , min{/c n , fc n }} (0 < n < M). 

Here 1 G Z M+1 is the vector whose entries are all 1. The minimum is to be 
interpreted entrywise, meaning that for x,y G R M+1 , 

min{x,y} = (min{x , yo}, min{xi, yi}, . . . , min{x M , Vm}) G R M+1 . 

Proof. The constraints ( 6.12[ ) imply that V(k) C Vq x • • • x Vm- It remains to 
show that any b G Vq X • • • X Vm belongs to V(k), which entails showing that 
there exists a scattering sequence p such that («(p), /3(p)) = (k,b). It suffices 
to construct a tree representing such a p, as follows. Given b G Vo x • • • x Vm, 
place k n vertices, consisting of b n branch points and k n — b n leaves (in any 
order), at depth z n , for < n < M, and place a root at z_i (considered 
as a branch point). Let iV denote the largest index such that k^ ^ 0. For 
< n < N, 

1 < &n-l < min{fc n _i, k n }, 

so there exists a surjection /„ from the k n vertices at z n onto the 6 n _i branch 
points at z n -\. The tree representing p is completed by drawing an edge 
from each vertex v at z n to f n {v). I 



6.3.1 A simple formula for the weight 

The tree representation facilitates deriving a simple formula for the weights. 
The following lemma uses multi-index notation, whereby given a vector S = 
(So, Si, ... , Sm) and an integer vector d = (do, . . . , dM), 

M 

s d =nst. 

n=0 
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Lemma 6.2 Let p £ Sj/ fe a scattering sequence in an M -layer model 
(t,R), and set (k,b) = (re(p), /3(p)). Then 

w (p) = (-Rf- b R k - b T 2b . 

Proof. Consider the tree representing p. Observe that each instance of 



w n = Rj in (6.6) corresponds to a unique leaf at Zj (see Figure 17). Since 



there are k n — b n leaves at z n this results in a total contribution of R k ~ b . 

Let v be a branch point at Zj having d v edges to vertices at Zj+i. Observe 
that precisely d v — 1 of these edges correspond to an instance of w n = —Rj 
in (6.6), and every occurrence of w n = —Rj arises this way. The sum total 
of numbers d v — 1 over branch points v at Zj is simply k n+ \ — b n = k n — b n , 
making for a total contribution over all depths of (—R) k ~ b . 

Finally, each instance of transmission from the jth layer to the (j + l)st 
layer in p corresponds to a vertex in the tree representing p which is not 
a leaf, i.e., to a branch point at Zj\ and, since the path p starts and ends 
at Z-\ every such transmission has a corresponding return transmission in 
the opposite direction, from the (j + l)st layer to the jth layer. There are 
bj branch points at Zj, and each of these corresponds to two transmissions 
across the boundary at Zj, making for a total contribution to w(p) of T 2b . 
Since every w n is covered by one of the above cases, the lemma follows. I 

6.4 An explicit formula for the impulse response 



A formula for the coefficients a(R,k) in (6.8), defined as the summation 



(6.9), is now within easy reach. Recall that the binomial coefficient ( J for 
a pair of non- negative integer vectors x,y £ Z +1 , with y < x, is to be 
interpreted as 



M 

X ' 



n 



(The inequality y < x means that x — y has non- negative entries.) 

Theorem 6.1 Let (r,R) be an M-layer model for some integer M > 1, let 
k € &M; an d se t u = min{l, k}. Then 

a(R,k)= £ (^(IZl)^^ 11 ^ 2 ^ ( 6 - 14 ) 

where V{k) denotes the set of b € Z M+1 such that u <b < minjfc, k}. 
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Proof. The total amplitude resulting from scattering sequences having a 
given transit count vector k G £m is 

a(R,k)= W (P)- 

PGSm such that 

By Proposition |6, 1| and Lemma |6.2| the above sum may be rearranged as 



a(R,k) = E W (P) 



b£V(k) p£Sm such that 
(«(p),/3(p))=(fc,6) 

= Yl (#{p e s ^ I ( k (p)' /3(p)) = {-Rf~ b R k - b T 2h . 

b£V(k) 

Thus all that is required is to count the number of scattering sequences 
having a given transit count vector and branch point count vector. Equiv- 
alently, it suffices to count the number of corresponding trees — which is 
straightforward. Consider first the arrangement of vertices in a tree for 
which (k(p),/3(p)) = (k,b). At each depth z n , there are k n vertices of which 
b n are branch points and k n — b n are leaves. There are ( 6 n ) ways of arranging 
these from left to right, making for a total of 

!)-n(t) ™ 

possible vertex arrangements. (There is only one way to place the root at 
Z—i, which may be ignored.) 

For each vertex arrangement there are various possible edge arrange- 
ments, as follows. Each of the k n +i = k n vertices at z n +i must be con- 
nected by an edge to one of the b n branch points at z n , respecting the 
vertex ordering (so that edges don't cross). This is equivalent to choosing 
a 5 n -part ordered partition of the integer k n . If k n > 1, there are (j^Zi) 
possible choices; and if k n = then b n = and there is 1 = (^ n ) (empty) 

arrangement. Letting iV denote the largest index for which k^ ^ 0, the 
total number of edge arrangements is 

k Hi \ I — r ( 1 



b — u, 

n=0 



n:: < 6i6 > 



12 June 2012 



Peter C. Gibson 



Reflections in layered media 



43 



Combining (6.15) and (6.16) yields a total tree count of 

#{ P gS m |(k(p) ! /3(p))= Q 
completing the proof. 



Thus the coefficients a(R,k) in the expansion (6.8) coincide with the 



amplitude polynomials a(R, k) of Definition 2.2 proving Theorem 2.1 

7 Conclusions 



For generic models, Theorems 4.1 and |4.2| express precisely the local linear 



algebraic character of the model-data correspondence, and the inverse al- 



gorithm, Algorithm 5.2 gives a practical method by which to recover the 
model from the data. Figures [2j [5] and [7] convey the geometric view of the 
underlying correspondence between travel time vectors r and arrival time 
vectors a. There are several remarks to be made about the implications of 
this general picture; each of the sections below considers a particular issue. 

7.1 Modifying the boundary depth Z-\ 

The physical setup in which source and receiver sit at depth Z—\, above 
the boundary zq of the layered half-space, may be easily modified to suit 
a particular application without altering the essential results or even the 
formulas. For example, consider the situation of source and receiver at depth 
z-i = zq. One has to decide whether a pulse must cross the boundary zq 



to be detected by the receiver. If not, then the formulas in Theorem 2.1 
remain valid by simply setting tq = and by replacing each occurrence of 
the transmission coefficient 



To = y l - Rl 

by 1. Similar modifications suffice to adapt the results to other physical sce- 
narios that might more accurately model a particular experimental modality, 
such as marine seismic etc. The results of the present paper are also rele- 
vant to non-seismic applications such as single source-receiver tussue sensing 
adaptive radar [22], [12"] . 



7.2 Non-generic models 

What happens as the normalized travel time vector f approaches the bound- 
ary of one of the cells in Figure [5] or its higher dimensional analogues? There 
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are two subcases to consider. If f crosses one of the orange lines — so that 
the corresponding lattice set remains constant — then the enumeration 
function permutes the values of two k,k' £ 27~ M . That is, there is an n 
such that on one side of the line ip{k) = n and ip(k') = n + 1, while on 
the other side ip{k') = n and ip(k) = n + 1, this being the only change 
in ip. On the line itself ip ceases to be injective, mapping k and k' to the 
same value. The second case is where f crosses a black line. In this case 
the lattice set S7 M gains or loses an element (so that the dimension of 
goes up or down accordingly). On the black line there is a transit count 
vector k ^ k M such that ip(k) = ip(k M ). The data corresponding to non- 
generic model associated with f on one of the lines need not determine a 
unique model, as in Theorem |3.1| But this does not mean that the inverse 
problem cannot be dealt with. (After all, a number of historical treatments 
are restricted to the case of equal travel times, a particular non-generic sce- 
nario.) If an extra quantity, i max = |t|, is included as part of the data, then 
the model can be recovered. However, there is a qualitative difference: the 
travel time inversion and amplitude inversion no longer fully decouple as 



they do in Alogorithm 5.2 Instead, the two must be interwoven, resulting 
in a slower algorithm that has additional steps. Thus the generic case is not 
only typical, it is also qualitatively simpler than the non-generic case. 

7.3 Coarsely layered media 

In working with finite precision, the meaning of generic must be appropri- 
ately adapted. Referring again to Figure[5j the normalized travel time vector 
f must be sufficiently far from the set of line segments to be considered as ef- 
fectively generic. This rules out the regions of the triangle close to the upper 
boundary where the lines are clustered ever more densely. Roughly speak- 
ing, there is a remaining part of the triangle that corresponds to effecively 
generic travel times — interior to the larger cells and sufficiently far from the 
boundary. A similar restricted set of reflection vectors R, sufficiently far 
from the algebraic hypersurfaces on which a(R, k) = (for appropriate k 



as in Section 4.1), completes the description of effectively generic models. 
It seems appropriate to refer to the layered media corresponding to such 
generic models as being coarsely layered. Not only is the distribution of 
depths restricted for such models, so is the number of layers. (This is be- 
cause the relative volume of cells in the higher dimensional analogues of 
Figure [5] goes down with increasing dimension.) Without attempting to 
quantify the notion more precisely, there is a regime of coarsely layered 
media in which the deterministic approach of the present paper constitues 
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a practical theory. For finely layered media where there is no effectively 
generic model, the stochastic methods of [8] are more appropriate. In any 
case the results presented here provide tools — both in the form of explicit 
formulas and, more broadly, a geometric perspective — with which to assess 
the utility and limitations of a deterministic approach to the 1-D acoustic 
reflection problem. 
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